Matrix- valued Boltzmann Equation for the 
Hubbard Chain 



Martin L.R. Furst 1,2 ' a , Christian B. Mendl 1,b , Herbert Spohn 1,3 ' 

August 15, 2012 

Abstract 

We study, both analytically and numerically, the Boltzmann transport 
equation for the Hubbard chain with nearest neighbor hopping and spa- 
tially homogeneous initial condition. The time-dependent Wigner func- 
tion is matrix-valued because of spin. The H-theorem holds. The nearest 
neighbor chain is integrable which, on the kinetic level, is reflected by 
infinitely many additional conservation laws and linked to the fact that 
there are also non-thermal stationary states. We characterize all station- 
ary solutions. Numerically, we observe an exponentially fast convergence 
to stationarity and investigate the convergence rate in dependence on the 
initial conditions. 
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1 Introduction 



The Hubbard model is a simplified description of interacting electrons moving 
in a periodic background potential, see [1, 2, 3] for introductory literature. We 
are interested in the dynamics of the Hubbard model in the regime of small 
interactions, which is conveniently described by kinetic theory, following the 
pioneering work of Pcicrls [4] , Nordheim [5] , and Uehling, Uhlenbeck [6] . From 
the point of view of kinetic theory the Hubbard model has unusual features. The 
non-interacting model has a doubly degenerate band, which - because of spin - 
makes the Wigner function 2x2 matrix-valued. In addition the hamiltonian is 
invariant under a global SU(2) spin rotations. On the kinetic level this property 
is reflected by an exceptionally large set of conserved quantities. We refer to [7] 
for a recent experimental realization through ultracold atoms in an optical lattice 
under conditions where also kinetic theory is applied. 

As one would expect, even the matrix- valued Boltzmann equation satisfies 
the H-theorem. The goal of our note is to achieve - beyond mere entropy in- 
crease - a quantitative and more detailed understanding of the approach to 
stationarity. The Boltzmann equation consists of the sum of an effective hamil- 
tonian plus a dissipative collision term, both with cubic nonlinearity. At such 
generality numerical simulation is not an easy task. Therefore we concentrate 
on the Hubbard chain with nearest neighbor hopping and on-site interaction. 
In addition we assume spatial homogeneity. Our simulations use 64 grid points 
in momentum space, which still allows for easy exploration. At this stage the 
reader might wonder, why on the kinetic level in one dimension there are any 
collisions at all. This will be explained in due course, as well as the difference be- 
tween the nearest neighbor integrable model and the non-integrable next nearest 
neighbor case. 

Let us start with the underlying hamiltonian and the resulting kinetic equa- 
tion. The electrons are described by a spin-^ Fermi field on Z with cre- 
ation/annihilation operators satisfying the anticommutation relations 

{a* a {x),a T (y)} = 8 xy 5 ar , {a a (x), a T (y)} = 0, a*(y)} = (1.1) 

for i,i/eZ, <t, r € {j, I}, and {A, B} = AB + BA. The hamiltonian reads 

H= a(x-y)a*(x)-a(y) + ~Y( a *W- a W) 2 - ( L2 ) 

Here a*(x) ■ a(x) — a|(x) Of (a;) + a|(x) a\,{x). a is the hopping amplitude, with 
the properties a{x) = a(x)*, a(x) = a(—x), and A is the strength of the on-site 
interaction. Our notation emphasizes the invariance under global spin rotations. 
For the Fourier transformation we use the convention 

f{k) = Yf{x)^ 2mk ' x - (1-3) 



Then the first Brillouin zone is the interval T — [— 1, 1] with periodic boundary 
conditions. The dispersion relation oj(k) = a(k) and, up to a constant, in 
Fourier space H can be written as 
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with k = k\ + k 2 — k 3 — k± mod 1 and d 4 fe = dki dk 2 dk 3 d&4. 

To arrive at the kinetic equation, we assume that the initial state of the 
chain is quasifree, gauge invariant, and invariant under spatial translations. It 
is thus completely characterized by the two-point function 

(al(k)a T (k')) = 6(k - k')W aT (k). (1.5) 

It will be convenient to think of W(k) as a 2 x 2 matrix for each k 6 T. Then, in 
general, W{k\)W{k 2 ) ^ W(k 2 )W(ki) and every argument of standard kinetic 
theory has to be reworked. By the Fermi property we have < W(k) < 1 as a 
matrix for each k. In particular, W can be written as 

W(k) = ]T e a (k)\k,a){k,a\, (1.6) 

where \k,a) for a G is a fc-dependent basis in spin space C 2 and e a are 

the eigenvalues with < e a < 1. 

At some later time t the state is still gauge and translation invariant, hence 
necessarily 

(al(k,t)a T (k',t)) = 6(k - k')W aT {k,t). (1.7) 

In general W(t) is a complicated object, but for small coupling A the quasi- 
free property persists over a time scale of order A -2 , a structure which allows 
one to obtain the kinetic equation by second order time- dependent perturbation 
theory. More details can be found, e.g., in [8, 9, 10]. Here we only write down 
the resulting Boltzmann equation 

~W(k,t) =C c [W](k,t) +C d [W](k,t) =C[W](k,t), (1.8) 
at 

which has the structure of an evolution equation and has to be supplemented 
with the initial data W(k, 0) = W(k). 
The first term is of Vlasov type, 

C c [W](k,t) := -i[H eS (k,t),W{k,t)}, (1.9) 

where the effective hamiltonian H c ^(k, t) is a 2 x 2 matrix which itself depends 
on W . More explicitly, 

#eff,i = dfc 2 dfc 3 dfc 4 S(k) V (i) 

x (W 3 Wi - W 2 W 3 - W 3 W 2 - tr[Wi]W 3 + tr[W 2 ]W 3 + W 2 ) . (1.10) 

Here and later on we use the shorthand W = 1 — W, W\ = W(ki,t), -ff e ff,i = 
H c ff(ki,t), uj_ = w(fci) +u{k 2 ) — co(k 3 ) — ui{k4). Since W is 2 x 2 matrix-valued, 
tr[-] is the trace in spin space. Finally V denotes the principal part. Since the 
k 3 , fei integration can be interchanged, H e e = H* s , as it should be. 

There are many different ways to write the collision term d- We choose a 
version which separates the various contributions into gain and loss term. Then 

C d [W]x = 7T ( dk2dk 3 dk i S(k)5(tj)(A[W} 1 234 + A[W}* l234 ), (1.11) 

JT3 
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where the index 1234 means that the matrix _4[W] depends on fci, fe, k 3 , and 
k<±. Explicitly 

A\W]um = ~W 4 W 3 W 2 + W A tr[WW 2 ] 

- {WiW 2 - WiW 3 - W 3 W 2 + W 4 tr[W 3 ] - Wi tv[W 2 ] + tv[W 2 W 3 }}W-i 

(1.12) 

with the first two summands the gain term and {...}Wi the loss term. The gain 
term is always positive definite, as implied by the inequality 

Atr[BC] + Ctx[BA] - ABC - CBA > (1.13) 

valid for arbitrary positive definite matrices A, B, C . Thus if an eigenvalue of 
W(k,t) happens to vanish, the gain term pushes it back to values > 0. A 
similar argument can be made for W{k,t) 1 implying the propagation of the 
Fermi property [10], to say: if at t = one has < W(k) < 1, then the solution 
to (1.8) also satisfies < W(k,t) < 1. 

In our contribution, we report on a numerical solution of the kinetic equa- 
tion (1.8), emphasizing the approach to stationarity. To provide an outline, 
in Sec. 2 we establish a few general properties of (1.8), (1.10), (1.11) . They 
hold for arbitrary lo and also for the obvious extension of (1.8) to d dimensions. 
In particular, we show that the entropy production a[W] = ^.S[W] has the 
property a > 0. The thermal state Wfb (the Fermi-Dirac distribution) satisfies 
C[W"fd] = an d hence also ct[Wfd] = 0. But to list all stationary solutions 
of (1.8) is not an easy task in general. 

In Sec. 3 we restrict ourselves to the Hubbard chain with nearest neighbor 
hopping, i.e., 

uj(k) = 1 -cos(27rfc). (1.14) 

The first task is to discuss the kinematically allowed collisions, in other words 
the solutions of w = together with k = mod 1. The nearest neighbor model 
has a special symmetry through which a large set of further stationary states, 
beyond the thermal ones, can be found. On the kinetic level, this reflects the 
integrability of the underlying quantum hamiltonian. In Sec. 4 our numerical 
procedure is outlined and in Sec. 5 it is used to study the dynamics for repre- 
sentative initial Wigner functions. 



2 General properties of the Hubbard kinetic equa- 
tion 

To emphasize generality, for this section only, we consider Z d as underlying 
lattice. Hence kj € T d with periodic boundary conditions. The SU(2) invariance 
of H is reflected by 

C[U* WU] = U*C[W]U (2.1) 

for all U G SU(2). Hence if W(k,t) is a solution to (1.8), so is U* W{k, t) U. 
Also hermiticity is propagated in time, i.e., if W(0) = W(0)* , then also W(t) = 
W(t)*, which follows from 

C[W}*=C[W*}. (2.2) 
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Furthermore the Fermi property, < W(t) < 1, is propagated in time, see [10] 
for details. 

There are two conservation laws 



• spin 

d 

dt J T 

• energy 



dkW(k,t)=0, (2.3) 



— / dkLj(k)tx[W(k,t)]=Q. (2.4) 



dt 



The proof uses the symmetrization of the integrand. One can interchange the 
variables ki O fc 2 , k% ■<-» fc 4 and also the pairs {fci,A; 2 } O \k^,ki}. For the 
energy, one then picks up the integrand lj and hence the factor wiS(w) = 0. 

The next general property is the H-theorem. Since |A| <§; 1, locally the state 
is free fermion. On the kinetic level, the entropy of the state W is then defined 
as 

S[W] = - I dki (tr[Wi log W x ] + tr[Wi log Wi]) . (2.5) 

Hence the entropy production is given by 
d f 

a[W] = —S[W] = - / dMr[(logVFi -logW^ClWh]. (2.6) 
dt J T d 

The H-theorem asserts that 

a[W] > for all W with < W < 1. (2.7) 
To establish (2.7), for each k we write 

W(k)= ^{k)Pa(k) (2.8) 

with eigenvalues < £ CT (fc) < 1 and orthogonal eigen-projections P a (k) — 
\k,a){k,a\ with (fc,cr|fc,</) = 8 acr i. As before, we use a shorthand as Pj = 
P^ikj), Ej = e aj [kj) and £ CT = £ ffl)0a|0Sf<u . Verting (2.8) into (2.6), one 
obtains 

a[W]=n / d 4 fc^(fc)5(u;) Y] (logei - logei)(£i£ 2 e3£4 - ei^^^) 

x (tr[PiP 3 ]tr[P 2 P 4 ] + tr[PiP 3 ]tr[P 2 P 4 ] - tr[P P 3 P 2 P 4 ] - tr[P 4 P 2 P 3 P 1 ]) 

= 7T / d 4 fc <5(ft)5(w) V" (£i£ 2 e 3 £ 4 - £i£ 2 £ 3 £ 4 ) log(e"i/ei) 
J(T*)* „ 

1 1 2 

X (fci,CTi|fc 3 ,CT 3 )(fc 2 ,Cr 2 |fc4,cr 4 ) - (fci , CTi | fc 4 , CT 4 ) (fc 2 , CT 2 |fc 3 , CT 3 ) • 



(2.9) 



We interchange 1 o 2, 3 O 4 and (1, 2) o (3, 4). Then 

£i£ 2 £ 3 £ 4 \ 
£ie 2 e 3 e 4 y 



a [ w ] = 7 / d4fc <H£)£fe>) (ei£2e 3 e4 - £i£ 2 £ 3 £4) log 

4 J (Td) 4 ^ 



1 1 2 

x (A;i,cri|fc 3 ,(j 3 )(A; 2 ,<T 2 |fc 4 ,a 4 ) - (fci,ai|fc 4 ,cr 4 )(fc 2 ,(T 2 |A; 3 ,cr 3 ) > 0, (2.10) 
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since (x — y) \og{x/y) > 0. 

Stationary states are defined by 

C[W] = 0, (2.11) 

which obviously implies cr[W] — 0. Physically one would expect thermal equi- 
librium to be included in the stationary states. On the kinetic level thermal 
equilibrium is defined by the Fermi-Dirac state 

WMk)= J2 (e^^-^ + l)' 1 ^)^!. (2-12) 
<re{t,4.} 

which is characterized by the inverse temperature /3, the two chemical potentials 
Hf, fj-l for the spin occupations, and some fc-independent spin basis \a), ft, jif, 
fil € K. Indeed, it is easily checked that C[Wfd] = 0. 

With this background information, one develops the following rough pic- 
ture on the approach to stationarity. The initial state determines a special, 
/c-independent basis \a) through 

dkW(k)= E *k)H (2-13) 

By (2.3) this basis is preserved in time. Thus it is natural to expand W(k,t) in 
this special basis. Approach to the thermal state would mean 

lim (a\W(k,t)\(j') = for cr^o-' (2.14) 

t— f oo 

and 

lim (a\W(k,t)\a) = (^M*)-"-) + for a e {t4}- (2.15) 

Since, by (2.3), the integral over the eigenvalue is conserved, one concludes that 

dk (V^W-aO + l) 1 , a e {t, i}- (2-16) 



dJfc(^"W-^ + l] 
Correspondingly, by (2.4), for the average energy 

e= / dkuj(k)tT[W{k)]= [ dk V oj(k) ffJXpW-i*) + A' 1 . (2.17) 

Both equations determine the parameters j3, [if, m from the initial W . One 
has < Ef, ei < 1 and min^ u>(k) < e < maxt to(k). Then the map (e, Ef, ej.) to 
(/3,/j,f,ni) is one-to-one. 

Implicitly our argument assumes that the set of stationary states equals the 
set of thermal states. But this might fail if there are not enough collisions, 
which could very well be the case in low dimensions. The issue of characterizing 
all stationary states has been accomplished only partially, see [11] for results 
towards this goal. On the other hand we still succeed in listing all stationary 
states and their domain of attraction. As to be discussed in the following section, 
for the Hubbard chain with nearest neighbor hopping the stationary states are 
not exhausted by the thermal ones. 
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3 Nearest-neighbor Hubbard chain 
3.1 Collisions 

w(k) 




Figure 1: The dispersion relation w(k) of (1.14). 



We return to the Hubbard chain with nearest neighbor hopping (1.14). Fig. 1 
visualizes w(k) for k £ [— |, |]. The first task is to investigate the kinematically 
allowed collisions defined by 6(k)8(w). The momentum conservation k = 
mod 1 allows us to eliminate one k- variable, say k 2 = k% + k^ — k\ mod 1. 
Inserted into energy conservation lo — and using some trigonometric identities, 
one arrives at 

lo = 4 sin(7r(fci — k 3 )) sin(7r(fci — k^)) cos(7r(/c3 + k^)). (3-1) 

Fig. 2 visualizes lu for fixed k\ = ||. From (3.1), we conclude that the collision 
manifold has a solution path k% + k^ = g (and thus also k\ + k 2 = | ) denoted 
7diag in Fig. 2, besides the "trivial" solutions k 3 — k\ (denoted 71) and k^ = ki 
(denoted 72). 

In what follows, we investigate the integral (1.11) of the dissipative collision 
operator along the paths 71, 72, and 7diag- Using the invariance of the 
integral (1.11) under k$ O ki, we may interchange W3 O W4. Then the 
integrand in (1.11) can be decomposed as 

A[W] 12U + A[W]* 12M = A,uad[W] 123 4 + A tr [W} 123 4 (3.2) 

with 

Ajuad[VF] 123 4 := -W1W3W2W4 - WiWiWsWx + W 1 W 3 W 2 Wi + WiWzWzWx, 

(3.3) 

Ar[lF]i234 := (WiWa + WWi)tr[WW 4 ] - {W X W 3 + W 3 Wi)tx[W 2 W4]. 

(3.4) 

Inspection of (3.3) immediately reveals that *4 qua d[lF]i22i = along 72 since 
(ki,k 2 ) = (/c 4 ,fc 3 ). Moreover, we also have ^4 qU ad[ll / ]i2i2 = along 71 with 
(ki,k 2 ) = (^3,^4), which can be checked by expanding (3.3). In other words, 
-4quad[lF]i234 contributes only along 7 diag . 
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Figure 2: Contour (blue straight lines) and gradient (green vectors) of the energy 
conservation uj = for fixed ki = If and after eliminating &2- The diagonal 
blue line 7di a g is precisely the contour k 3 + = ^. The vertical and horizontal 
blue lines, 71 and 72, are the contours k^ = fci and £4 = k\, respectively, pi 
marks the intersection of ji with 7diag for i = 1,2. 



The situation is different for the term *4tr[W^]i234: while Ar[W / ]i2i2 = 
along 71 by direct inspection of (3.4), it is (in general) non-zero along 72 and also 
along 7diag- In summary, for evaluating the dissipative collision integral (1.11) 
we have to integrate >Ur[W] along 72 and both AniadP^] and ^Ur^] along 

7diag- 

As a side remark, the solution path 7di ag is special for the nearest neighbor 
dispersion relation (1.14). If we add to (1-14) a small next-nearest neighbor 
term, then 71 and 72 persist and 7di ag gets somewhat deformed. In addition, a 
new collision channel opens up, as illustrated in Fig. 3 for the dispersion relation 

^nnn(fc) = w(fc) — ^ cos(47rfc) = 1 — cos(27t/c) — ^ cos(47r/c). (3.5) 



3.2 Stationary solutions 

The collision paths 71, 72, and 7diag have special symmetries, from which one 
can guess the form of stationary solutions beyond the thermal one. They have 
the same structure as the Fermi-Dirac state, but with w(fc) replaced by a more 
general function /. One finds 

W st (k)= M*)k)H M*) = (e /W - a -+l)~\ (3.6) 

<re{t4} 
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Figure 3: Contour (blue straight lines) and gradient (green vectors) of the energy 
conservation for a next-nearest neighbor model with w nnn (A:) of (3.5) and fixed 

k, - 23 

«1 - 64- 



where / is a real-valued, 1-periodic function satisfying f(k) = — /(§ — k), a a € R, 
and \a) is an orthogonal basis, independent of k. 

As discussed in the Appendix, (3.6) characterizes the entire set of stationary 
solutions. The next step is to identify the domain of attraction for W st , in 
other word to study the map from the initial W to W s t ■ Here we can follow the 
strategy described at the end of Sec. 2. 

We first note that there are many energy-like quantities which are conserved. 
Let g : T -> R with g(k) = -g{\ - k). Then 

~^dkg(k)tT[W(k)]=Q, (3.7) 

which generalizes the energy conservation (2.4). (3.7) again follows by an ap- 
propriate interchange of the integration variables ki, ■ ■ ■ , k^. 

By substituting g(k) = 5(k — k') — 5{k — | + k') for arbitrary k' € T, one 
concludes that 

h(k) = tr[W(fc)] - tr[W(| - k)] (3.8) 

is pointwise constant for each k € 1. Assuming that the initial W converges to 
a stationary state of the form (3.6), it must hold that 

h ( k ) = E ((e /(fc) ^+l) _1 -(e- /(fc) - a -+l) _1 ). (3.9) 
ae{t,J.} 

Equivalently, as in Sec. 2, the spin conservation law requires that the eigen- 
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values e a in (2.13) are equal to 

dfc (e /(fc) - a - + l) 1 . (3.10) 

We claim that (3.9) and (3.10) uniquely determine / and a ai or more specifically, 
that the map between 

tx[W{k)]-tt[W{\-k)l\k\<\, 0<£ t , £i <l (3.11) 

and 

f(k) with f{k) = -/(§ - k), \k\ < \, at, a x (3.12) 

is one-to-one. In particular, to a given W one can associate a unique W st of the 
form (3.6). 

Proof. By a short calculation, (3.9) can be written as 

h(k) = sinh(/(fc)) ( coshat+ \ oshm + coshai+ \ oshm ) ( 3 - 13 ) 
and (3.10) as 

sinh a, 



i \ cosh a a + cosh f(k) ^ 



with interval of integration I := [— j, j}. We define a generalized "free energy 
through 

H{f,a^,a i ) = J^dk ^ log ( cosh a CT + cosh /(fc)) . 



(3.15) 



?e{t,l} 

The map (/, a^,a^) i— > _ff is strictly convex. Furthermore 



d v [ai sinha CT 1 

- — H = / dfc — — — — — — s (j ~ (3.16) 

oa a J l cosh a a + cosh/ (fc) 2 

and 

sinh/(fc) _ 

*/(*0 ~ ^ cosha a + cosh/(A;) ~ W ' 1 j 

Thus the map from above can be viewed as Legendre transform from the first 
set (3.11) to the second set of variables (3.12). Since H is convex, the map is 
one-to-one. □ 



4 Numerical Procedure 

4.1 Mollifying the collision operators 

Dissipative collision operator. We have to make sure that 5(uj)S(k) is a well- 
defined prescription. For this purpose we eliminate ki and, using (3.1), obtain 




dfc 4 S(uj) = \d k Mk^ kl r 1 = (2tt |sin(27rfc 3 ) - sin^fcx)!) 1 . (4.1) 
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Likewise along 7di ag it holds that 



&k A 5(uj) = |5 fe4 w| fe4=1/2 _ fcE 



I 1 = (2tt |sin(27rfc 3 ) - sin(27rfci)|) 1 



(4.2) 



Idi; 



Considering the subsequent integration over k% in (1.11), the critical point ky, = 
\ — k\ (marked p 2 in Fig. 2) would lead to infinities in general. (Integrating along 
Tdiag across the point p\ (k^ = fci) is possible since -4[VF]i234 +-4[W]i234 is zero 
at that point, as explained above). As mollification we choose the substitution 



(2?r |sin(27rfc 3 ) - sin(27rfci)|) 1 (4tt 2 ( sin(27rfc 3 ) - sin(27rfci)) 2 + e 



-1/2 



(4.3) 

with some finite e > 0. Concretely, we use e = | for the simulations. 

Conservative collision operator. The integral (1.10) for the conservative collision 
operator C c differs from dissipative integral (1.11), since there is only a single 
delta distribution S(k). Thus, we can eliminate fc 2 — fc 3 + k^ — k± as for the 
dissipative case, but still have to integrate over both fc 3 and k^. 




(a) 



(b) 



Figure 4: (a) The term 1/uj as function of k% and k4, for fixed fci = 23/64. (b) 
The "mollified" version w/(w 2 + e 2 ) with e = \ is free of singularities. 



The integral (1.10) is defined as Cauchy principal value with respect to 
1/uj. Fig. 4a illustrates this term in dependence of &3 and k<± (compare also 
with Fig. 2). While the Cauchy principal value exists for continuous W(k), the 
numeric calculation is rather demanding and we resort to a mollifying procedure 
as for the dissipative collision operator. Concretely, we substitute 



1 

u 



(4.4) 



with finite e > (in our case e = |). Fig. 4b shows the right-hand side, in direct 
comparison with the unmollified version. Note that C c could be defined via the 
integral (1.10) with the replacement (4.4), and then letting e — » 0. 
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4.2 Solving the Boltzmann equation 

In order to solve the Boltzmann equation (1.8) numerically, we discretize the k 
variable by a uniform grid 

k, = ] -, j = 0,...,n-l (4.5) 

with n = 64 in our case. We have chosen the interval [0, 1] instead of (equiva- 
lently) [—1/2, 1/2] simply for convenience. Note that due to periodicity, W(l, t) = 
W(0,t), so the point k = 1 is not required. We use the trapezoidal rule to ap- 
proximate the integrals (1.11) and (1.10) of the dissipative and conservative col- 
lision operators, respectively. Note that this approach is particularly suited for 
analytic period functions. Moreover, considering the 2-dimensional integral of 
the conservative collision operator, we ensure that the variable ki = k% + k± — k\ 
mod 1 is a grid point whenever k\, k^ and k^ are grid points, in distinction from 
other integration rules with non- uniform points. 

We solve the Boltzmann differential equation (1.8) for the time variable by a 
Strang splitting (or symmetric Trotter splitting) technique: denoting the (fixed) 
timestep by At, we combine an explicit midpoint rule step for the dissipative 
part with the time evolution operator for the conservative part: 

X(kj,t) = e- iJr -*<*i'*> At/2 W(kj,t) e ar ««(*i'*) A */ 2 , j = 0, . . . , n - 1, (4.6) 



Y(kj,t) = X{kj,i) + At C d 



X(t) + ^C d [X(t)] 



(kj), j = 0, ...,n-l, 

(4.7) 

W(kj,t + At) = e-^tf^'*) At / 2 Y(kj,t) e^CM) j = , . . . , n - 1, 

(4.8) 

where H' eff depends on Y(t). The midpoint rule has order 2, while the time 
evolution operator eT lH At / 2 (-)e lH At / 2 has only order 1. Thus, the complete 
integration scheme has order 1. As advantage, the time evolution operator 
preserves matrix symmetry. For the simulations, we use At = 1/16, and the 
overall simulation time interval runs from t = to varying upper limit t = 
15,..., 45. 



4.3 Cost analysis 

Considering a single time step, the most expensive part is the evaluation of 
the conservative collision operator C c in (4.6) and (4.8), i.e., the 2-dimensional 
integral (1.10) after eliminating /c2- (The dissipative collision operator Cd re- 
quires only a one-dimensional integration.) For the uniform discretization with 
n points in each direction, this scales like 0(n 2 ). One time step requires the 
evaluation of this integral for n different k\ points, thus the overall cost is 0(n 3 ). 

On a Intel Core i7-740QM Processor (6M cache, 1.73 GHz) without using 
parallelization, one time step takes approximately 90 s (Mathematica 8 imple- 
mentation), so a complete simulation is approximately 6h. Note that the per- 
formance could be easily increased by a C/C++ implementation and making 
use of parallelization. 
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5 Simulation results 



High-temperature state. Fig. 5 illustrates a high-temperature Fermi-Dirac equi- 
librium state Wfd (2.12) where /3 = 1CT 4 , ^ = 10 4 and m = -10 4 . We have 



W FD (k) 
1.0 r 



0.8 



0.6 
0.4 



0.2 



-0.4 



-0.2 



0.0 



0.2 



0.4 



Figure 5: Diagonal matrix entries of a high-temperature equilibrium state Wfd- 
The state is (almost) independent of k. 

chosen the initial state (see Fig. 6) by W(k, 0) = Wfd(&) + V(k) with V(k) a 



W(k,0) 

1.0 r 



eig(W(k,0)) 
1.0 r 





(a) matrix entries of initial W(k, 0) 



-0.4 -0.2 0.0 0.2 0.4 
(b) eigenvalues of W(k, 0) 



Figure 6: (a) Initial state W(k, 0) = W FU (k) + V(k) with V(k) defined in (5.1). 
The blue and green curves show the real diagonal entries, and the red and 
magenta curves the real and imaginary parts of the off-diagonal |t)(4-| entry, 
respectively, (b) Corresponding eigenvalues of W(fc,0) in the interval [0, 1]. 



rotation of the Pauli a z matrix and subtracting the constant matrix r/18: 



V(k) = ^e- 2ntTk cr z e 2 ' TiTk - ^, 



1 



T = (T. r — <7„ 



V(k) satisfies 



/ dkV(k) = and tr[V(fc)] = 0, 
Jt 



(5.1) 



(5.2) 



for all k G T such that W(0) matches Wfd in terms of the spin and energy 
conservation laws (2.3) and (2.4). 
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Fig. 7 illustrates the convergence to the equilibrium state Wfd- Interest- 
ingly, we observe that the off-diagonal entries converge slower than the diagonal 
entries, but an analytic explanation of this effect is still lacking. 



l|W(t)-W FD || 




S[W FD ]-S[W(t)] 
0.1 



decay rate: 
-0.337939 



decay rate: 
-0.648043 




2 4 6 8 10 12 14 

(a) convergence in Hilbcrt-Schmidt norm 

||diag(W(t)-W FD )|| 
0.1 
0.001 



2 4 6 8 10 12 14 

(b) entropy convergence 



|(W(t)-W FD ) 12 | 
1 r 




decay rate: 
-0.337867 



2 4 6 8 10 12 14 

(c) convergence of the diagonal entries 



2 4 6 8 10 12 14 

(d) convergence of the off-diagonal entries 



Figure 7: Convergence of the initial W(0) (Fig. 6) to the high-temperature 
equilibrium state VFfd (Fig. 5) as semi- logarithmic plot (blue). The decay rate 
is the slope of the fitted red dotted line. 

Figure 8 displays the Bloch vectors f(k,t) € K 3 of W(k,t) parametrized by 
k, i.e., 

W{k,t) = ~{l+r>(k,t)-a), <j=(<j x ,a y ,a z ). (5.3) 

The dark blue curve shows the initial r(k, 0) and lighter blue curve shows r(fc, |). 
As time progresses, the initial curve straps to almost a single point, since Wpi)(k) 
is almost independent of k. 

Low-temperature state. Fig. 9 illustrates a low-temperature Fermi-Dirac 
equilibrium state Wfd (2-12) with /3 = 7, m = j| and //j, = y|. In this case, 
for given Wpu(k) the variational freedom for the initial W(k, 0) with the same 
symmetries as M / fd(^) is strongly restricted. Similar to the high-temperature 
state, we define W(k,0) = W FB {k) + V(k) (see Fig. 10) with 



V(k) = \ ( 
Again, V(k) satisfies 



64sin(7r(A;-3/4)) 2 _ e ~64 sin(7r(fc-l/4))' 



— 2iri a x k pP-fti cr x k 



/ dkV(k) = and tt[V(k)] = 0, 
Jt 



(5.4) 



(5.5) 
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(a) (b) 

Figure 8: Bloch sphere representation (dark blue curve) of the initial W(k, 0) 
(Fig. 6), parametrized by k and viewed from 2 perspectives. The light blue 
curve shows the corresponding Bloch curve of W(k,t) for t — 1/2. Finally, the 
curve for the high-temperature equilibrium state Wpuik) is indiscernible from 
a single point at the tip of the z-axis arrow. 

W FD (k) 







// °- 8 




// 06 




// °- 4 




// °- 2 





-0.4 -0.2 0.2 0.4 

Figure 9: Diagonal matrix entries of a low-temperature equilibrium state Wfd 
(/3 = 7, fPf = jg, m = 

for all k £ T. We observe that the convergence to the equilibrium state (Fig. 11) 
is slower than for the high-temperature state in the previous paragraph. (Note 
that the simulation time interval is now [0,45] as compared to [0, 15].) 

Degenerate chemical potentials. We consider a Fermi-Dirac equilibrium state 
with degenerate chemical potentials /i^ = as illustrated in Fig. 12. As initial 
state W(k,0), we set W(k,0) = W FB {k) + V{k) (see Fig. 13) with V(k) taken 
from (5.1). As illustrated in Fig. 14, there is no indication that the convergence 
changes due to the degeneracy. Two time snapshots of the eigenvalues of W(k, t) 
are shown in Fig. 15. They have a peculiar shape, and converge to the diagonal 
entries of Wfd, as expected. 

Non-thermal stationary state. For this example, we start from an (rather 
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eig(W(k,0)) 







1 °' 8 




// °' 6 




// 0,4 




// 02 





(a) initial W(k, 0) 



-0.4 -0.2 0.2 0.4 

(b) eigenvalues of W(k, 0) 



Figure 10: (a) Initial state W(k, 0) = W FD (k) + V(k) with V(k) denned in (5.4). 
The blue and green curves show the real diagonal entries, and the red and 
magenta curves the purely imaginary off-diagonal entries, (b) Eigenvalues of 
W(k, 0) in the interval [0, 1]. 



l|W(t)-w FD || 
0.050 



S[W FD ]-S[W(t)] 



decay rate: 
-0.0501206 





decay rate: 
-0.10088 



0.001 



10 20 30 40 

(a) convergence in Hilbert-Schmidt norm 

||diag(W(t)-W FD )|| 



10 20 30 40 

(b) entropy convergence 



|(W(t)-W FD ) 12 | 
0.1 




10 20 30 40 

(c) convergence of the diagonal entries 



decay rate: 
-0.227469 



10 20 30 40 

(d) convergence of the off-diagonal entries 



Figure 11: Convergence to the low-temperature equilibrium state Wfd (Fig. 9) 
as semi-logarithmic plot (blue). The decay rate is the slope of the fitted red 
dotted line. For this example, we observe that the off-diagonal entries converge 
much faster than the diagonal ones. 



arbitrary) initial 
W(k,0) = 



-l e -cos(47r(fe- 7 )) + 1 I sm ( e 2™fc) 

5 V \ sm(e- 27Tik ) |erf(cos(27rfc)) + f + arctan(sin(27rfc - |)) + 4 

16 " (5.6) 




Figure 12: Diagonal matrix entries of an equilibrium state Wfd with /3 = 1 and 
same chemical potentials /x-j- = //j, = 1. 



eig(W(k,0)) 




k 



(a) initial W(fc,0) (b) eigenvalues of W(fc,0) 

Figure 13: (a) Initial state W(k, 0) = Wfu (k) + V{k) with Wpuik) proportional 
to the identity matrix (see Fig. 12), and V(k) defined in (5.1). The blue and 
green curves show the real diagonal entries, and the red and magenta curves 
the purely imaginary off-diagonal entries, (b) The eigenvalues of W(k,Q) are 
non-degenerate, different from the equilibrium state Wpuik). 



illustrated in Fig. 16 (where 7 is the Euler gamma constant), and then deter- 
mine the stationary, non-thermal state W st {k), via the /-function described in 
Sec. 3.2. Fig. 17 illustrates both / and W s t- Next, we run the numerical sim- 
ulation of the time evolution, which should converge to the predicted W s t(k). 
Fig. 18 indeed verifies the convergence to W st (k). 

For special cases, we have checked that the asymptotic decay rate is almost 
independent of the initial conditions. This strongly suggests that the collision 
operator linearized at W s t has a spectral gap. 
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l|W(t)-W FD || 

1 




decay rate: 
-0.350924 



0.001 



2 4 6 8 10 12 14 
(a) convergence in Hilbert-Schmidt norm 



S[W FD ]-S[W(t)] 




decay rate: 
-0.690033 



2 4 6 8 10 12 14 

(b) entropy convergence 



Figure 14: Convergence to the equilibrium state Wfd with degenerate eigenval- 
ues (Fig. 12) as semi-logarithmic plot (blue). The decay rate is the slope of the 
fitted red dotted line. 



eig(W(k,t)), t=l/2 eig(W(k,t)), t=2 




1 k 1 

-0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 

(a) eigenvalues at t = i (b) eigenvalues at t = 2 



Figure 15: Two snapshots showing the convergence of the eigenvalues to the 
equilibrium state Wfd (red, same as Fig. 12) with fif = [i± = 1. 



eig(W(k,0)) 



W(k,0) 


1.0 


0.8 




0.8 


0.6 




f \ 0.6 






\ / 0.4 


o.r 







-Q>^ t 0.2 
(a) initial W(k, 0) 



-0.4 -0.2 0.0 0.2 0.4 
(b) eigenvalues of W(k, 0) 



Figure 16: (a) Initial state W(k,0) defined in (5.6). The blue and green curves 
show the real diagonal entries, and the red and magenta curves the real and 
imaginary parts of the off-diagonal |t)(4| entry, respectively, (b) Eigenvalues of 
W(k, 0) in the interval [0, 1]. 
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-0.4 -0.2 0.0 0.2 0.4 

(b) 



Figure 17: (a) The /-function (blue) calculated from tr[W(k,0) - W[\ - fc,0)] 
(dashed) and the fitted "chemical potentials" = —0.617485 and — 
0.0578622. The initial W(k,0) is defined in (5.6). (b) Resulting stationary 
state W st (k) (3.6) given by / and a^, a^. 



S[W sl ]-S[W(t)] 




2 4 6 8 10 12 14 



(a) entropy convergence 



HW(t)-W st || 

1.000 

0.500 




2 4 6 8 10 12 14 



(b) convergence in Hilbert-Schmidt norm 



Figure 18: Convergence to the calculated W st (Fig. 17) as semi-logarithmic plot 
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6 Conclusions 

The kinetic equation for the Hubbard model, in general, has two hardly inves- 
tigated features (i) the Wigner function is 2 x 2 matrix- valued, (ii) the micro- 
scopic SU(2) invariance implies additional conservation laws. We investigated 
here the chain with nearest neighbor hopping, which is an integrable model, [1]. 
The Boltzmann transport equation reflects integrability by an infinite number 
of conserved quantities and non-thermal stationary states. We established the 
H-theorem and classified all stationary states. Adding a next-nearest neigh- 
bor coupling seems to destroy all conservation laws beyond spin and energy 
which indicates that now the stationary solutions are exhausted by the thermal 
Fermi-Dirac Wigner functions. 

In the spatially homogeneous case we observed numerically an exponentially 
fast convergence to the predicted stationary state, both for the diagonal and 
off-diagonal matrix elements with roughly comparable decay rates. The decay 
at low temperatures is slower than at high temperatures, as one would have 
expected. In principle, asymptotic decay rates can be computed from the lin- 
earized collision operator. 

Physically of great interest would be to better understand the spatially in- 
homogeneous situation. For example one could imagine to have in each half of 
the chain a thermal state with the same temperature, but with different spin 
orientations. In principle, this could be handled by kinetic theory. One only 
would have to add in the kinetic equation the transport term w'(fc) d/d x . Nu- 
merically, such a problem is more demanding than the one studied here, but, 
at least in one dimension, still in reach. Another challenging problem would 
be to study energy transport through the chain. Our results point towards the 
validity of Fourier's law. 

A Characterization of stationary solutions 

Proposition 1. Let a[W] be as defined in (2.6). If < W < 1, then the 
solutions to zero entropy production, 



are necessarily of the form (3.6). 

Remark. As noted by J. Lukkarinen, further zero entropy and stationary so- 
lutions are obtained by setting one eigenvalue of W identically = 0, 1, and the 
other eigenvalue arbitrary. 

Proof. On the one hand, if W is of the form (3.6), then cr[W] = follows 
by inserting. On the other hand, let a[W] = 0. We set fc = (^1,^2,^3,^4), 
er = (<7i, (72, 03, 04), d k = dfci d&2 dk 3 d&4 and define 



with the £i as in Sec. 2, furthermore 

1 1 2 

G(k,a) = \(k 1 ,o- 1 \k 37 a 3 )(k2, 0-2^4,0-4) ~ (fci, crilh, cr 4 )(fc 2 , cr 2 \k 3 , cr 3 ) . (A.3) 



a[W] = 0, 



(A.l) 




(A.2) 
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Then 

a[W} = j[ d 4 W(fc)^)Vf(fc !( r)G(fe,(T) (A.4) 

according to (2.10). Since all terms are non-negative, 

F{k,a)G{k,a) =0 (A.5) 

must hold for all u and all k E 72 U 7diag (see Fig. 2). On 71 one has F = and 
no extra information can be extracted. 

F has the structure (x — y) log(-), which is zero only iff x = y, equivalently 

iff 

los (St) = 108 (!) + loe (!) - be (!) - los (!) = °- < A - 6 » 

Defining the collision invariants as 

*' (t, = loe (!i)- <AJ) 

condition (A. 6) reads 

*<r a (fci) + $ ff2 (fa) - $ ff3 (fc 3 ) + $ CT4 (fc 4 ). (A.8) 

Note that the labeling of eigenvalues £f(k), £i(k) and corresponding eigen- 
vectors is arbitrary. Thus w.l.o.g. we can assume that 

(Ai,tl*2,t) 7^0 and thus (h,i \k 2 ,i) ^ (A.9) 
for all k 1 ,k 2 ET. 

Consider the contour 72 {k\ — £4, k 2 = £3) for cr =^|tl- In this case, the 
second term on the right side of (A. 3) vanishes, and thus 

G(k, nil) = I (fci , t I k 2 , t) (h , I I k 2 , 1) 1 2 > (A. 10) 

by construction (A.9). Therefore (A.5) forces F(fc, fiti) = on 72. Equa- 
tion (A.8) becomes after rearranging terms 

$ t (fci) - $4.(fci) - $ t (fc 2 ) - ^(k 2 ). (A.ll) 

Since variables are separated, both sides of (A.ll) must be constant, i.e., 

$ t (fc)-^(fc) = c (A.12) 

for a fixed eel and all k € T. 

Next, we establish that the basis \k,a) has to be ^-independent up to a k- 
dependent phase, which can be chosen such that \k, cr) = \a) with \\), \\) a fixed 
basis in C 2 . If c = in (A.12), then $ T = $ i; and it follows that W{k) = e(k) 1. 
In particular, one can set \k, a) = \a). In the other case, c ^ 0, consider the 
contour 72 for cr =nil : 

F(fe,tm) = $ T (fci) + $ T (fc 2 ) - - ^(fci) - 2c ^ 0, (A.13) 
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where we have used (A. 12) for the second equality. Thus (A. 5) requires that 
G(fc,tt!l) = on 72. Inserted into the definition (A. 3) yields 

(k 1 ,t\k 2 ,l)(k 2 ,t\k 1 ,l)=0 (A.14) 

for all k\,k2 € T. Since the vectors |fc,t) and \k,\) are an orthonormal basis of 
C 2 for each fixed k, (A.14) is equivalent to 

(k 1 ,t\k 2 ,l) =0 (A.15) 

for all k\,k 2 G T. Keeping k 2 fixed, this means that |fei,t) = const up to a 
phase, and similarly \ki,\) = const. W.l.o.g. the phase factor can be set to 1, 
leaving invariant the projectors P a {k) — \k,a){k,u\. In summary, \k,o~) = |cr) 
and G(jk,tr) = G(tr). 

As final step, consider 7di ag for er =^m. By direct inspection G(titl) = 1) 
thus (A. 5) requires F(k, t-j-tl) = 0. (A. 8) for k 2 = \ - ki and fc 4 = \ - k 3 
becomes 

$ t (fci) + $ 4 (| - fci) = $ t (fc 3 ) + - k 3 ). (A.16) 
Since variables are separated, both sides must be constant, i.e., 

$ t (fc) + - jfc) = const (A.17) 
for all fe G T. Combined with (A. 12), we obtain 

$ CT (fc) + - fc) = const (A.18) 
for cr =f, One concludes that $ a is necessarily of the form 

**(*:) = - *a with f a {k) = -f a {\-k) (A.19) 

for all A: € T and some a CT € K. Plugging into (A. 12), one deduces that f^(k) — 
fl(k) — const and, since f a (\) — 0, it follows that f\{k) = fi{k) = f(k) 
independent of a. Summarizing, we arrive at 

$ CT (fc) =f{k)-a a . (A.20) 

Solving (A.20) and (A. 7) for e a (k) leads to the claimed form (3.6). □ 

Corollary 2. Under the constraint < W < 1, all stationary solutions, i.e., 
all solutions to C[W] = 0, are precisely of the form (3.6), 

W st (k)= K(k)\a)(a\, X a (k) = (e^' + l)" 1 (A.21) 

^e{T,;} 

with f(k) = -/(| - k) for all k e T. 

Proof. Each W s t of the form (A.21) satisfies C[W] — 0, which can be checked 
by inserting W st into C[W]: specifically, the commutator (1.9) defining C c [W s t] 
vanishes since H c g and W s % are diagonal. The dissipative collision operator 
Cd[Wst] is zero due to the symmetry properties of 7di a g and the fact that f(k) = 
~f{\ ~ k ). On the other hand, let W be a solution of C[W] = 0. Then 

~W(k,t)=C[W](k,t)=0, (A.22) 

and, in particular, 

o-[W] = ^-S\W] = 0. (A.23) 
at 

According to Proposition 1, W is of the form (A.21). □ 
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